GlacierUpdate Subroutine

SUBROUTINE GlacierUpdate   & 
  !
  ( time, mask, rainfall )       

IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT (IN) :: time
TYPE (grid_integer), INTENT (IN) :: mask

!Arguments with intent(inout):
TYPE (grid_real), INTENT(INOUT) :: rainfall  !rainfall (liquid precipitation) (m/s) 

!local declarations:
INTEGER (KIND = short) :: i, j, is, js
REAL (KIND = float) :: meltRate
REAL (KIND = float) :: tAir
REAL (KIND = float) :: cmelt !!snow melt coefficient (m/s/°C)
REAL (KIND = float) :: tmelt
REAL (KIND = float) :: icewe_tminus1 !!icewe at previous time step
REAL (KIND = float) :: ddx	!!actual cell length (m)
REAL (KIND = float) :: ds !!thickness of the saturated depth (m)
REAL (KIND = float) :: cellWidth !!cell width (m)
REAL (KIND = float) :: ijSlope  !!local slope (m/m)
REAL (KIND = float) :: qin, qout  !!input and output water discharge in snowpack
REAL (KIND = float) :: area !!cell area (m2)
INTEGER (KIND = short) :: currentDOY !!current day of year
!------------------------------end of declarations-----------------------------  

!check if parameter update is required
CALL GlacierParameterUpdate (time)

!snow ice transformation
currentDOY = DayOfYear (time, 'noleap')
IF ( dtSnow > 0 .AND. currentDOY == doySnowTransform ) THEN
    DO i = 1, mask % idim
        DO j = 1, mask % jdim
            IF (mask % mat(i,j) /= mask % nodata ) THEN
                icewe % mat (i,j) = icewe % mat (i,j) + swe % mat (i,j)
                swe % mat (i,j) = 0. 
                waterInIce % mat (i,j) = waterInIce % mat (i,j) + waterInSnow % mat (i,j)
                waterInSnow % mat (i,j) = 0.
            END IF
        END DO
    END DO
    
END IF

!initialize melt grid
iceMelt = 0.


!compute inter-cell lateral water flux
QinIce = 0.
QoutIce = 0.
DO i = 1, mask % idim
    DO j = 1, mask % jdim
        IF ( mask % mat (i,j) == mask % nodata ) CYCLE
        
        !set downstream cell  (is,js)
        CALL DownstreamCell ( i, j, flowDirection % mat(i,j), is, js, ddx, &
                            flowDirection ) 
        
        IF ( CheckOutlet (i, j, is, js, flowDirection) ) CYCLE
      
        !thickness of the saturated depth
        ds = waterInIce % mat (i,j)
        
        !cell width
        cellWidth = ( CellArea (mask, i, j) ) ** 0.5
        
        !local slope
        ijSlope = ( dem % mat (i,j) - dem % mat (is,js) ) / ddx
        
        IF ( ijSlope <= 0. ) THEN
            ijSlope = 0.0001
        END IF
        
        QoutIce % mat (i,j) = cellWidth * ds * iceConductivity % mat (i,j) *  ijSlope
      
        !output becomes input of downstream cell
        QinIce % mat (is,js) = QinIce % mat (is,js) + QoutIce % mat (i,j)
        
    END DO
END DO

!loop across cells
DO i = 1, mask % idim
    DO j = 1, mask % jdim
        
        !skip nodata
        IF ( mask % mat (i,j) == mask % nodata ) THEN
            CYCLE
        END IF
        
        !compute potential melt rate
        tAir  = temperatureAir % mat (i,j)
        cmelt = meltCoefficientIce % mat (i,j) / 1000. / 86400.
        tmelt = meltTemperature % mat (i,j)
        SELECT CASE ( meltModel % mat (i,j) )
        CASE (1) !degree-day temperature based
            meltRate = DegreeDay ( tAir,  tmelt,  cmelt )
        CASE DEFAULT
            CALL Catch ('error', 'Glacier', 'melt model not implemented: ', &
                         argument = ToString (meltModel % mat (i,j) ) )
        END SELECT
        
        ! actual melt rate, limited by available ice 
        ! in the previous time step
        meltRate = MIN ( meltRate, icewe % mat (i,j) / dtIce )
    
        !update icewe
        IF (dtSnow > 0 ) THEN
            IF ( swe % mat (i,j) <= 0.) THEN !no protection by snow
                icewe_tminus1 = icewe % mat (i,j)
                icewe % mat (i,j) = icewe % mat (i,j) - meltRate *  dtIce
                IF ( icewe % mat (i,j) < 0. ) THEN !meltrate too high
                    icewe % mat (i,j) = 0.
                    iceMelt % mat (i,j) = icewe_tminus1
                ELSE
                    iceMelt % mat (i,j) = meltRate *  dtIce
                END IF
            ELSE  ! snow protects ice from ablation
                iceMelt % mat (i,j) = 0.
            END IF
        ELSE  !no protection by snow
            icewe_tminus1 = icewe % mat (i,j)
            icewe % mat (i,j) = icewe % mat (i,j) - meltRate *  dtIce
            IF ( icewe % mat (i,j) < 0. ) THEN !meltrate too high
                icewe % mat (i,j) = 0.
                iceMelt % mat (i,j) = icewe_tminus1
            ELSE
                iceMelt % mat (i,j) = meltRate *  dtIce
            END IF
        END IF
        
        !update water in snow
        area = CellArea (mask, i, j) 
        qin = QinIce % mat (i,j) / area
        qout = QoutIce % mat (i,j) / area
        
        !lateral water flux occurs even when ice is covered by snow
        IF ( icewe % mat (i,j) > 0. ) THEN !rainfall contributes to water in ice
            waterInIce % mat (i,j) = waterInIce % mat (i,j) + &
                                        iceMelt % mat (i,j) +  &
                                         rainfall % mat (i,j) * dtIce + &
                                        ( qin - qout ) * dtIce
             rainfall % mat (i,j) = 0.
        ELSE 
            waterInIce % mat (i,j) = waterInIce % mat (i,j) + &
                                        iceMelt % mat (i,j) +  &
                                        ( qin - qout ) * dtIce
        END IF

        
        IF ( waterInIce % mat (i,j) < 0. ) THEN
            waterInIce % mat (i,j) = 0.
        END IF
        
    END DO
END DO

!write point icewe
IF ( time == timePointExport ) THEN
   CALL GlacierPointExport ( time )
END IF

RETURN
END SUBROUTINE GlacierUpdate